Sampling diffusive transition paths 
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We address the problem of sampling double-ended diffusive paths. The ensemble of paths is 
expressed using a symmetric version of the Onsager-Machlup formula, which only requires evaluation 
of the force field and which, upon direct time discretization, gives rise to a symmetric integrator 
that is accurate to second order. Efficiently sampling this ensemble requires avoiding the well-known 
stiffness problem associated with sampling infinitesimal Brownian increments of the path, as well 
as a different type of stiffness associated with sampling the coarse features of long paths. The fine- 
feature sampling stiffness is eliminated with the use of the fast sampling algorithm (FSA), and the 
coarse- feature sampling stiffness is avoided by introducing the sliding and sampling (S&S) algorithm. 
A key feature of the S&S algorithm is that it enables massively parallel computers to sample diffusive 
trajectories that are long in time. We use the algorithm to sample the transition path ensemble for 
the structural interconversion of the 38-atom Lennard- Jones cluster at low temperature. 
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I. INTRODUCTION 

Transition path sampling (TPS) is a successful ap- 
proach to characterizing rare event processes for which 
a priori knowledge of the location of, or even existence 
of, a single well-defined transition state is not available. 
It addresses the rare event problem by focusing compu- 
tational resources on the small subset of trajectory space 
that exhibits dynamically interesting events. The de- 
velopment of the approach into a computationally effec- 
tive tool is due largely to the efforts of Chandler and 
coworkers ) 1 ' 2 ' 3 ' 4 ' 5 ' 6 and it is now widely used to study 
the mechanisms and rates of reactive processes. 

The TPS method prescribes Metropolis Monte Carlo 
sampling of the transition path ensemble, which includes 
only those dynamical trajectories of a given length in 
time that begin in a prescribed "reactant" region and ter- 
minate in prescribed "product" region. A central feature 
of the TPS method is the preservation of the Boltzmann 
distribution through detailed balance. 

As in all Metropolis algorithms, TPS involves a pro- 
posal step, in which a new dynamical trajectory is con- 
structed, and then an acceptance/rejection step. New 
trajectories can be proposed in a variety of ways. One 
approach is to employ open-ended trajectory algorithms, 
which generate new trajectories by directly integrat- 
ing the dynamics of the system. Examples of this ap- 
proach include the widely used "shooting" and "shifting" 
techniques i 2 - An alternative approach employs double- 
ended trajectory algorithms, in which the dynamical tra- 
jectories are treated as extended geometrical objects, 
or chains-of-stateSji£i£ and sampled using path integral 
Monte Carlo methods. 

In practice, two factors have contributed to the pre- 
dominance of TPS applications that employ opened- 
ended trajectories Firstly, for processes that react 



by way of a single barrier-crossing event, the open-ended 
shooting and shifting algorithms are very effective in pro- 
viding decorrelated members of the transition path en- 
semble, thus rendering a double-ended approach unneces- 
sary. Secondly, the path integral sampling problem that 
arises in a double-ended formulation of TPS has appeared 
to be prohibitively difficult. 

However, strong motivation for the development of 
double-ended TPS algorithms does exist. Many complex 
systems, including many biologically relevant systems, 
undergo transition processes that involve slow, diffusive 
barrier crossing dynamics and exhibit the possibility of 
metastable intermediates. These conditions are known 
to cause difficulty for the open-ended TPS algorithms, 5 
but they are naturally overcome in double-ended meth- 
ods which enable the trajectory endpoint positions to be 
constrained to the reactant and product regions. Fur- 
thermore, recent developments in path integral Monte 
Carlo methods and the availability of massively parallel 
computers suggests that a double-ended formulation of 
TPS can be made tractable. 

In the present work, we consider the issue of sam- 
pling double-ended trajectories for diffusive dynamics. 
Going beyond the coarse-graining efforts of Elber and 
coworkers we seek the development of an efficient, 
ergodic sampler of the transition path ensemble. We 
show in Section II that the transition path ensemble for 
diffusive dynamics can be expressed using a symmetric 
Onsager-Machlup formula. In Section III, we show that 
discretization of this formula gives rise to a symmetric nu- 
merical integrator that requires only first derivatives of 
the potential energy function and exhibits second order 
accuracy. We draw contrast with the standard Onsager- 
Machlup formulajii the discretization of which gives rise 
to the stochastic Euler integrator that is only first or- 
der in accuracy. In Section IV, we address the path in- 
tegral Monte Carlo sampling problem that accompanies 
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a double-ended TPS approach. We begin by establish- 
ing the need for fast sampling^ the fine features of the 
stochastic process. We also argue that additional meth- 
ods are needed to efficiently sample the coarse features 
of trajectories that are long in time. We then introduce 
the sliding and sampling algorithm that alleviates this 
coarse-feature sampling problem and gives rise to a nat- 
ural parallelization scheme. In Section V, we utilize the 
S&S algorithm to sample the ensemble of reactive trajec- 
tories for the low-temperature interconversion of the 38- 
atom Lennard- Jones cluster between its octahedral and 
icosahedral minima. 



II. THE SYMMETRIC ONSAGER-MACHLUP 
FORMULA 

The time-evolution of relatively large particles in a 
highly viscous thermal bath is described by overdamped 
stochastic dynamics. In one dimension, these dynamics, 
which are also known as Brownian or diffusive dynamics, 
are governed by the law 



xt+j^V'ixt) = W t . 



(1) 



Here, Xt and it the particle position and velocity at time 
t, and V{xt) is the external potential energy function. 
The quantity 7 is the friction coefficient, which is related 
to the diffusion coefficient (D = l/((3j)) via the recip- 
rocal temperature (3. The white noise Wt is formally 
expressed as the time derivative of a Wiener process, or 
Brownian motion, of variance 2Dt. Extension of Eq. fl}, 
and the rest of this presentation, to higher dimensions is 
trivial. 

The object of primary interest in this study is the dis- 
tribution of transition paths,— £ 



P AB (x;t) = I A (x )P(x-;t)I B (x n ), 



(2) 



where P(x; t) is the equilibrium distribution of trajec- 
tories that obey Eq. ([1]) and pass through points x = 
{xq, Xi, . . . , x n } in subsequent timesteps of At = t/n. 
Ia{x) is the indicator function for the subset of config- 
uration space associate with the reactant species, and 
Ib(x) is the corresponding product indicator function. 
From the Markov property of diffusive dynamics, it fol- 
lows that in the canonical ensemble 



P(x;t) 



J| K(x k ,x k _i; At), 



(3) 



fc=i 



where K(x',x;t) is the conditional probability that the 
particle arrives at point x' , given that it started at point 
x some time t earlier. 

K(x',x;t) obeys the Smoluchowski equation, 



d d 
° K(x',x;t)=D^ 



dt 



_d 

dx 



7 + (3V'(x') 



K{x',x;t)\, 



which is the partial differential equation that corresponds 
to Eq. ((T|). However, by reweighting the conditional prob- 
ability, 

G(x', x- 1) = e + ^ v ^- v ^ 2 K(x', x- 1), (5) 

Eq. ((4|) is cast in the self adjoint form of the Bloch equa- 
tion, 



d_ 

dt 



G(x',x;t) 



D 



d 2 



dx' 2 



V eS (x') 



G(x',x;t), (6) 



with 



V eS (x) = D 



t V > [x f _ t V »( x) 



(7) 



The solution for the conditional probability is then im- 
mediately known from the Feynman-Kac solution for the 
Green's function of the Bloch equation^ 



G{x\x-t)^ D fe^[~^f o 



V'{W u fdu 



D/3 



V"{W u )du\, (8) 



with the symbol E^ 2 ^*'* meaning an average over all 
Brownian paths of variance 2Dt that start at x and end 
up at x' in time t. 

In practice, it will be useful to eliminate the evaluation 
of second derivatives of the potential function (Eq. ([8])). 
We thus perform a final manipulation. Using that 

Km (W tk+1 -W tk ) 2 -^^ 2Dt \W tk+1 ~W tk f =2DAt k , 

we consider the difference between forward and backward 
Ito integrals, obtaining 



/ V'{W U ) ■ d»W u - / V'(W U ) ■ df\V u = hm V [V'(W tk+1 ) - V'(W tk )] (W tk+1 - W tk ) 
Jo Jo At ^% =0 

= lwa Q ^2V"(W tk )(W tk+1 - W tk ) 2 - 2Dma Q J2V"(W t JAt k =2D V"(W u )du. (9) 

i — n 1 — n **Q 
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Substituting the last equality in the Eq. (J5J) , we finally 
obtain the symmetric Onsager-Machlup formula 

G(x', x- 1) = E^<) exp {-^ J V'{W u fdu 

+ ^J*V'(W u )-(d b -df)W u y (10) 

The symmetric Onsager-Machlup formula in Eq. ()10D . 
which to our knowledge is an original contribution of this 
work, is so called because of its clear relationship to the 
standard Onsager-Machlup formula, 11 

K(x',x;t) = E^ ) exp(-^ fv\W u ) 2 du 



^J*V'(W u )-df\Y u y (11) 



However, unlike Eq. (|11[) . Eq. (|10[) is symmetric in x and 
x' . (That Eq. (|10p is symmetric is obvious from the in- 
spection of Eq. ((9| ; that Eq. (jTTJ) is not symmetric follows 
from its dependence on a forward Ito integral and the lack 
of equality of forward and backward Ito integrals.) Al- 
though the two solutions are formally equivalent, it will 
be shown in the next section that the symmetric formula 
gives rise to more stable numerical implementation. 

Before proceeding, we note that the Feynman-Kac 
equation is also symmetric in x and x' and promises good 
numerical stability, but we will currently avoid evaluat- 
ing second derivatives of the potential energy function 
and reserve exploration of this avenue for a later date. 



III. SECOND ORDER DISCRETIZATION OF 
THE SYMMETRIC ONSAGER-MACHLUP 
FORMULA 

Discretization of the standard Onsager-Machlup for- 
mula gives rise to Euler's integrator, iA a first-order in- 
tegration scheme. In this section, we will show that 
discretization of the symmetric Onsager-Machlup for- 
mula gives rise to a more stable, second-order integration 
scheme. 

A symmetric discretization can be obtained from the 
symmetric Onsager-Machlup formula by utilizing the fi- 
nite difference appearing in Eq. ^ for the difference of 
Ito integrals and the trapezoidal quadrature technique 
for the Riemann integral. Thus, letting tk = kt/n for 
all < k < n, Wk = t/n for all 1 < k < n — 1, and 
wo = w n = t/(2n), we obtain the approximation 



G(x,x';t)*E^ 



+ j E^Wfc) " V '( w t k -Mw tk - WVJ \ . (12) 



exp ■ 



(3 2 D 



J2w k V'(W t , 



k=0 



Since the joint distribution of the Wiener process is 
known explicitly, the right-hand side can be reduced to a 
finite-dimensional integral as follows. Define the quantity 



p a 2(x,x') 



1 



-(x'-x) 2 /(2 a 2 ) 



V2 



net* 



(13) 



and let 



a 2 = 2Dt. (14) 
Then, setting xq = x and x n = x' , 

G(x,x';t)ta I dxi--- I dx n - 1 p a 2 /n (x,x 1 )p l7 2 /n (x 1 ,X2) 



■Pa 2 /n{ x n-l,x')ex.p ■ 



(3 2 D 



J2^kV'{x k ) 2 (15) 



k=0 



+ t y^}V'( x k) ~ V'(xk-l)](xk ~ X k -l) 



k=l 



Finally, the right-hand side of Eq. (I15|) can be ex- 
pressed in a way that is more familiar to those working 
with imaginary-time path integrals. G(x,x';t) has the 
form of a Lie- Trotter product 



G{x,x']t) ps / dxi ■ ■ ■ / dx n -\G{)(x,xi;t/n) 
Jr Jr 

xG (xx,x 2 ;t/ri) . . .GQ{x n - 1 ,x';t/n) (16) 

obtained from the "short-time" approximation (remem- 
ber that a 2 = 2Dt) 

&o(x,x ;t) =p a 2{x,x )exp<^ 

+ ^[V'(x')-V'(x)}(x'-x)Y (17) 

Eqs. ©, ©, ©, and $T7J) fully specify the discretized 
distribution of transition paths employed in this study, 



I A (x Q )e-^^ 



Y]_ G Q (xk,x k -i;Atk) 



k=l 



lB{x n )e 2 



-§V{x n ) 



(18) 

In considering the symmetric formula, we have been 
motivated not only by the need to ensure detailed bal- 
ance, but also by the observation of Suzuki that the or- 
der of convergence of a short-time approximation of the 



form e 



-tA^-tB 



(which is generally one) is improved when 



the symmetric form e tA /' 2 e tB e tA/2 is utilized (to the 



k=l 



-tA/2 e -tB e -tA/2 

value of two, which is also the maximal value for such 
products; a summary of Suzuki's findings can be found 
in Ref. 14). A more general theory, which is also valid 
outside the domain of applicability of the Trotter the- 
orem, has been designed by one of us?^ That theory 
can be utilized both to verify the convergence order of a 
short-time approximation as well as to design improved 
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onesJ^ It also allows one to design optimal propagation 
techniques useful for grid implementations^ 

We shall utilize such a grid technique to verify numer- 
ically that the convergence order of the short-time ap- 
proximation given by Eq. (JTTJ) is two. The grid technique 
is a sparse version of numerical matrix multiplication. 
The choice of grid and propagation technique depends 
on the order of convergence that is desired. For example, 
if the convergence order of a spatially-continuous short- 
time approximation is two, then a grid technique of order 
two preserves the order of convergence, whereas a grid 
technique of order greater or equal to three preserves not 
only the order but also the convergence constant. We 
shall therefore utilize the fourth-order grid from Ref . [l6| 
to verify the convergence order. Once the order is ver- 
ified and found not to exceed two, one can utilize the 
second order grid technique for actual computations on 
a grid. The lower the order of the grid is, the better 
the sparsity and the speed of the propagation scheme. 
Grid techniques are generally limited to low dimensional 
systems. 

The physical system considered is taken from Ref. [lj. It 
is a two-dimensional two-channel problem described by 
the potential 

V(x, y) = (4(1 -x 2 - y 2 f + 2(x 2 - 2) 2 + {{x + yf - if 

+ ((x-y) 2 -l) 2 -2)/6. (19) 

In reduced units, the friction coefficient is 7 = 3, the 
inverse temperature is (3 = 8, whereas the total propaga- 
tion time is set to t = 60, twice the value utilized in the 
original reference. 

We shall verify the ability of the short-time approx- 
imation to preserve the Boltzmann distribution. Using 
Eq. (fT8")l . if we let Ia( x ) = Ib( x ) — 1 an d integrate over 
all coordinates, we should ideally obtain the value 

Z(fi) = f e- fw{x) dx. 

This is the case for the exact Green's function G(x, x';t), 
which leaves the distribution e~^ v ^/ 2 invariant. If the 
short-time approximation Go(x,x';t) is utilized instead, 
then the result after propagation, denoted here by Z n {j3), 
is normally different from Z{0). It converges to Z(f3) in 
the limit n — * 00. 

We want to establish how fast this convergence is. An 
easy way is to verify that the difference n 2 [Z n ((3) — Z((3)] 
converges to a constant different from zero (the conver- 
gence constant). As this is apparent from Fig. [TJ second- 
order convergence is numerically established. Because of 
the long times utilized, the number of path variables nec- 
essary for an accurate treatment of rare events is large. 
For our example, a number of slices equal to n = 2048 
produces a relative accuracy of —2.3% for the partition 
function. 
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FIG. 1: The instantaneous convergence constant for the par- 
tition function as a function of the number of time slices. 
The guiding line is the convergence constant determined for 
n = 16,000. 



IV. PATH INTEGRAL MONTE CARLO 
SAMPLING 

The distribution of trajectories defined by Eq. (fl8|) 
poses a path integral Monte Carlo (PIMC) sampling 
problem. We can thus employ existing PIMC sampling 
techniques as a starting point for sampling the ensemble 
of diffusive transition paths . 12 ' 17 However, the goal of 
sampling trajectories that are long in time gives rise to 
challenges that have not been fully addressed by current 
methods. These issues, and the techniques that we pro- 
pose for overcoming them, are discussed in the current 
section. 



A. Sampling Fine Features of the Path: 
The Fast Sampling Algorithm 

Define ly to be the characteristic lengthscale of the 
potential energy function, namely the lengthscale below 
which the forces arising from the potential function are 
slowly varying. Furthermore, define Ty = ly/2D, the 
typical time required for a free particle to diffuse a dis- 
tance ly . If we then consider a segment of length r from 
any diffusive trajectory such that r <§; ry, we see im- 
mediately from Eq. (|17jl that the distribution for this 
path segment is governed almost entirely by the Gaus- 
sian (free-particle) term. It follows that diffusive trajec- 
tories experience a nearly drift-free Brownian motion on 
short timescales, which implies that such trajectories are 
composed of nearly independent infinitesimal increments. 

As is well known , 12 ' 17 this near-independence of the 
infinitesimal increments must be respected in a sampling 
algorithm if the usual path integral stiffness problem is to 
be avoided. For example, sampling the path distribution 
with individual time-slice moves gives rise to an artifi- 
cial correlation between the increments (often blamed on 
the "harmonic springs" in the free particle term) which 



causes poor computational scaling with the number of 
time-slices nJZ 

A solution to this problem is to use collective Monte 
Carlo moves for the time-slices, or multi-slice moves. 17 
A general strategy for doing so is to express the diffusive 
trajectories using a random series representation, namely 
an infinite series of specialized functions, 18 

oo 

x + (x' — x)u 

i=l 

where < u < 1. The conditions for these specialized 
functions are dictated by the Ito-Nisio theorem.— Per- 
missible functions include the well known Wiener- Fourier 
(sine function) serieS ) 18 i 20 ' 21 as well as the Levy-Ciesielski 
series^ 2 - that will be discussed shortly. 

Random series representations have several important 
properties. In particular, the Aj(ti) are constructed such 
that their amplitudes diminish with increasing index i; 
the tail of the random series thus becomes nearly decou- 
pled from the potential function and describes the fine 
features of the path. Furthermore, in the special case of 
V'(x) = const, Eq. (|20|) yields the correct distribution of 
trajectories upon letting the eti be a set of independent 
identically distributed standard normal variables. To- 
gether, these properties ensure that the path variables in 
the random series representation become nearly decou- 
pled in describing the fine features of the path. That is, 
direct Monte Carlo sampling of the generates multi- 
slice moves that avoid introducing artificial correlation 
between the infinitesimal increments of the path. 

The Levy-Ciesielski random series (LC) representation 
is especially amenable to efficient path sampling. In par- 
ticular, the recently developed Fast Sampling Algorithm 
(FSA) utilizes the LC representation to provide, for a 
fixed trajectory time, arbitrarily high accuracy in sam- 
pling the fine features of the path at a computational scal- 
ing of only n log 2 n in the number of path variables^ 2 - We 
briefly review the Levy-Ciesielski random series represen- 
tation and the fast sampling algorithm, which are em- 
ployed in the numerical examples throughout this study, 
and we refer the reader to Ref. [l2j for a more complete 
discussion. 

In the LC representation, the path is expanded in the 
basis of Schauder functions, 

F kt3 {u) = 2-( fc - 1 )/ 2 F M (2 fc - 1 U - j + 1), (21) 

where k > 1, 1 < j < 2 fe ~ 1 , and 

( u uE (0,1/2] 
F 1A (u) = i 1 -u u e (1/2,1) (22) 
I otherwise. 

These functions, which are illustrated in Fig. 2, are or- 
ganized into fc-indexed layers that describe increasingly 
coarse features of the path. Each layer contains a set of 
functions which are only non-zero on the disjoint, 
open intervals (uk,j-i, Uk,j)i where 1 < j < 2 k ~ 1 and 
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FIG. 2: Plot of renormalized Schauder functions for the first 
three layers. 

uu.j = j2~( k ~ x \ Because of this property, we know that 
only one j-indexed function within a layer contributes to 
the expansion, such that 

2 k-i 

X/ bjFk,j{u>) = ^[2 fc - 1 ii]+i-P 1 fc,[2 fc - 1 «]+i(w) (23) 

for any sequence of numbers bi, b 2l • ■ ■ , b 2 k~i. (Here, [x] 
denotes the largest integer smaller than or equal to x, 
and for u = 1, the quantities 6 2 fc - 1 +i an d Fk.2 k - 1 +i{^) 
are defined to be zero.) We thus have the following LC 
representation of the path, 

00 

x+(x'-x)u+V2m^2a kt[ 2h-i u]+1 F k>[2 h-i u]+ i(u). (24) 
fe=i 

Using the LC representation, we sample the path 
variables a k j according to the distribution defined in 
Eq. (fT8|) . In doing so, we employ the Fast Sampling 
Algorithm (FSA), which prescribes that each path vari- 
able be sequentially updated using the Metropolis Monte 
Carlo algorithmic The FSA exploits two key properties 
of the LC representation to enables the efficient sampling 
of large numbers of time-slices. 

Suppose that the path is described using n + 1 time- 
slices (including endpoints), where n — 2 k , and con- 
sider a particular LC path variable akj such that k > 1. 
As is clear from the preceding paragraph, changing this 
variable involves moving only n k — 1 time-slices, where 
rik = 2 fc ~ k + 1 . A beautiful consequence of this fact is 
that, for path integrals of the Lie- Trotter product form, 
the change in the path integral action accompanying a 
change in ak.j can be obtained at only a fraction (« 2 1 ~ fe ) 
of the computational cost of evaluating the action for the 
total path. (This is an unusual feature for a random se- 
ries representation. For example, changing a single path 
variable in the sine function series requires reevaluating 
the action for the entire path.) We thus arrive at the 
key properties of the LC representation: (1) In sampling 
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a Lie- Trotter product form, the path variables that be- 
long to the same layer are strictly independent. (2) The 
cost of computing the changes in path integral action 
for a sequential sweep through all n path variables in 
the LC representation is only logarithmically greater (a 
factor k' times) the cost of doing so in the time-slice rep- 
resentation. Since it is proven that optimal statistical 
statistical efficiency of independent variables is obtained 
by a sequential Metropolis Monte Carlo scheme^ it is 
clear that the FSA exploits the first of these properties. 
The FSA also exploits the second property, which ensures 
that path variables can be introduced into the tail of the 
random series at a computational cost that scales as only 
n log 2 7i, to ideally solve the problem of sampling the fine 
features of the path. 

What is so special about the symmetric Onsager- 
Machlup formula that requires utilization of the FSA? 
The main source of concern is the appearance of Ito in- 
tegrals in the expression for the action. The difference of 
backward and forward Ito integrals appearing in Eq. (1101) 
would vanish if the path were a smooth function. In 
fact, if 4>(u) is piecewise smooth, then both the forward 
and the backward Riemann sums define a same Riemann- 
Stieltjes integral, that is, 



/ V'(W ut )d f cf>(u) = / V'{W ut )d b 4>(u) 
o Jo 

= [ V'(W ut )j)(u)du. (25) 
Jo 



Hence, the Ito integral in Eq. (fT0|) would vanish instead 
of converging to a quantity involving the Laplacian of the 
potential, as shown by Eq. ©. It is only the very fine 
details (the so-called noise) that prevents the difference 
of Ito integrals from vanishing. 

If, for example, a Lennard- Jones cluster at low tem- 
perature starts to suddenly exhibit unexplained liquid- 
like behavior during a path-integral simulation, the re- 
searcher should recall that there are two main reasons 
that the difference of Ito integrals vanishes accidentally: 
i) inadequate sampling of the fine details and ii) a num- 
ber of time slices that is too small for the chosen transit 
time (not enough fine details). If the difference of back- 
ward and forward Ito integrals vanishes or is too small, 
then the distribution of paths is controlled only by the 
square-norm of the gradient. The distribution tends to 
accumulate around the stationary points of the potential, 
where the gradient vanishes. However, the simulation has 
no way to discern whether these stationary points are lo- 
cal minima, saddle points, or local maxima. It is worth 
mentioning that the explicit inclusion of the Laplacian 
term through a numerically better-behaved Riemann in- 
tegral makes the Feynman-Kac formula less dependent 
upon the quality of sampling. 



B. Two Kinds of Sampling Stiffness 

A central observation of the current paper is the follow- 
ing: Although random series representations avoid stiff- 
ness in sampling the fine features of the path, they intro- 
duce stiffness in sampling its coarse features. By using 
multi-slice moves, random series representations intro- 
duce correlations between distant segments of the path. 
That is, changing a single path variable in a random 
series representation generally involves moving a large 
number of time-slices along the entire length of the path 
(imagine, for example, changing a coefficient in the sine 
function expansion). This is not a problem in sampling 
the fine features of the path, because the amplitude of 
the motion is small in comparison to ly - too small to 
"feel" the potential function. But the introduced cor- 
relation is a severe problem in sampling the coarse fea- 
tures of long paths. If the trajectories are long in time 
{t ^> Ty), then well-separated segments of the path will 
in fact be almost entirely uncorrelated - and should be 
sampled accordingly. By using a representation that ar- 
tificially correlates the well-separated segments - and at- 
tempts to move them over distances that are comparable 
to ly - we introduce stiffness. 

Stiffness in sampling the coarse features of the path 
is an obvious problem in the context of reactive trajec- 
tories, which travel from a reactant region to a distinct 
product region - and thus involve motion on lengthscales 
that are long compared to ly. However, we note that it is 
an inevitable problem for any path integral Monte Carlo 
application that utilizes a random series representation 
and involves long times (including imaginary times) . Re- 
gardless of the context, we would expect to encounter 
stiffness whenever dealing with potentials that have a 
characteristic lengthscale that is short in comparison to 
the lengthscale for free particle motion ((2Z?t) 1 / 2 in the 
current context, (flti 2 /m) 1 / 2 for imaginary time path in- 
tegrals). 

In the following example, we illustrate how stiffness 
associated with sampling the coarse features of the path 
appears in practice. We return to the two-dimensional, 
two-channel problem considered in the previous section. 
Using reduced units, we set the friction coefficient 7 = 3, 
the inverse temperature /3 = 8, and the trajectory time 
t = 60. We describe the path using n = 2048 time-slices, 
which corresponds to employing k' = 11 layers in the 
LC representation. The endpoints of the trajectory are 
held fixed at (x — ±1, j/ = 0). (Of course, additionally 
sampling the endpoints is trivial, but it is also superflu- 
ous when a sufficiently long t is employed.) Employing 
the FSA, the LC path variables are updated sequentially 
according to the Metropolis algorithm. We propose new 
path variables according to 

a fcj = a °k ,j + Aa fci3 -£, (26) 

where £ has probability density 2~ 1 (1 + a; 2 ) -3 / 2 . Ran- 
dom numbers having this distribution can be generated 
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TABLE I: Proposal distribution widths for the Metropolis 
Monte Carlo moves in the two-dimensional example.™ 
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a k is the layer index for the LC representation of the path, 
rik — 1 is the number of time-slices moved by changing the 
path variable 0,,/s, Aa*,^ is the equilibrated path variable 

distribution width, and Aa; is the corresponding distance in 
coordinate space. 

via the identity £ = (C — 0.5)/(C(l — C)) _1/2 , with ( uni- 
form on [0, 1]. The widths of the proposal distributions 
Acifcj are restricted to be the same for all path variables 
within a given fc-indexed layer, and the numerical values 
for these widths are tuned in an equilibration phase to 
yield approximately 40% acceptance. 

The results for the equilibrated proposal distribution 
widths Aflitj- are presented in Table |U We also pro- 
vide the corresponding distances in coordinates space 
Ax = (L>i) 1/2 2- fe / 2 Aa fc j, obtained using Eqs. fH) 
and {M)) , 

For layers fc = 9 — 11, the Aa^j reach a plateau, in- 
dicating that the tail of the random series expansion has 
decoupled from the potential energy surface and that the 
simulation is converged with respect to the number of 
path variables. In this tail of the expansion, we see from 
the Ax values that subsequent layers sample increasingly 
fine features of the path. 

However, for k < 9, the Aakj decrease with decreasing 
k, eventually falling off exponentially for the most coarse 
layers. This behavior indicates that path variables are 
no longer sampled from their free particle distributions 
and increasingly feel the presence of the potential energy 
surface. The onset of this decrease in proposal widths at 
k — 8 is reasonable, given that the characteristic distance 
traveled by a free particle in = 16 time-slices is about 
0.2 - a distance that is commensurate with the finest 
details of the potential surface. 

For layers k — 5 — 8 in Table [l] we note that although 
the distribution widths Aa^,j decrease, the correspond- 
ing distances moved by the time-slices remains large and 
relatively constant. This feature arises because LC layers 
with lower k indices correspond to larger lengthscales in 
configuration space. Despite the decrease in the proposal 
distribution widths, we see that the LC representation is 
still efficiently sampling the time-slices associated with 



these middle layers. 

However, for k < 5 we see an exponential decrease 
in both the proposal distribution widths and the corre- 
sponding time-slice moves in configuration space. This 
follows from the attempt of ever-coarser moves in the 
a random series representation and inevitably leads to 
coarse-feature sampling stiffness. 

C. Sampling Coarse Features of the Path: 
The Sliding and Sampling Algorithm 

In the preceding subsection, we recognized that ran- 
dom series representations, by utilizing multi-slice moves, 
introduce stiffness in the Monte Carlo sampling of the 
coarse features of long paths. We now offer a simple 
solution to this problem. We propose an algorithm in 
which the total path is divided into shorter fragments. 
Keeping its endpoints fixed, each fragment is represented 
with a random series and independently sampled using 
Metropolis Monte Carlo. The Metropolis-evolved frag- 
ment configurations arc then used to reconstruct the total 
path - such that all time-slices other than the fragment 
endpoints have generally moved. Finally, to ensure er- 
godicity, the total path is re-divided into a different set 
of fragments and the process is repeated. 

How does this approach avoid both coarse-feature and 
fine-feature sampling stiffness? Clearly, by using a ran- 
dom series representation for the fragments, we avoid the 
fine-feature sampling stiffness. Moreover, by fragmenting 
the path, we achieve ergodic sampling without employ- 
ing the multi-slice Monte Carlo moves that suffer most 
extremely from coarse-feature stiffness. 

This algorithm essentially generates a rejectionless dif- 
fusive dynamics for the fragment endpoints, which in 
turn define the domains for independent, Metropolis sam- 
pling. Since diffusive dynamics is itself just a sampling 
device, we are, in a sense, sampling the sampler. 

Our specific implementation of this idea utilizes the 
LC representation. Presume that the total path of time 
t is uniformly discretized into n + 1 time-slices, 

x= (a;o,ari,...,x n _ 1 ,x„). (27) 

We fragment the path into / segments containing rif + 1 
time-slices (including endpoints), where rif = 2 kf and 
frif < n. Each fragment can then be represented using 
an LC series with kf layers. Upon requiring that the 
fragments form consecutive segments of the total path 
(i.e. share endpoints), there are n+ 1 — fn t possible ways 
to fragment the path. For example, if n — 10, / = 2, and 
nf =4, the three possibilities are 

!(x ,X 1 ,X 2 ,X 3 ,X4), (x4,X 5 ,X 6 ,X 7 ,X S ), [xg , Xg , X W ] 
[xo,Xi] , (x 1 ,x 2 ,x 3 ,X4,x 5 ), (x 5 ,x 6 ,x 7 ,x$,x 9 ), [x 9 ,a;io] 
[X ,X 1: X 2 } , (x 2 ,X 3 ,X4,X 5 ,X(i), (x 6 ,x 7 ,x 8 ,x 9 ,x 10 ), 

where the segments included in parentheses are the frag- 
ments, and those included in square brackets are the left- 
over ends of the trajectory. When the path is re-divided 
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FIG. 3: Schematic illustration of a sliding move. Using 
Schauder tents, the LC representation for the path is shown 
in solid lines before sliding and in dotted lines after sliding. 
The segment of the path that corresponds to a particular frag- 
ment, before and after sliding, is shown using the solid and 
dotted boxes, respectively. 



between fragment sampling steps, we randomly choose 
from among these various possibilities, effectively slid- 
ing the consecutive fragment sampling domains along the 
stationary path. Each fragment sampling step involves a 
single sweep through the fragment LC path variables us- 
ing the FSA. We call this procedure sliding and sampling 
(S&S). 

The sliding move is illustrated schematically in Fig. [3] 
The top of the figure shows a segment of the total path, 
which extends out of view in both directions, as a func- 
tion of time. At the bottom, in solid lines, we see the LC 
representation for a particular fragmentation of this path, 
using lengths of n/ = 8. The LC layers are illustrated 
in the usual way with a nested structure of "Schauder 
tents." The solid box indicates the segment of the path 
that corresponds to a particular fragment. In dotted 
lines, the figure shows how the path is re-fragmented in 
a sliding move. Although the path is unmoved, its repre- 
sentation in the LC path variables changes as we slide the 
fragment endpoints to the right by one time-slice. The 
dotted box illustrates how the domain over which the 
path is sampled by the corresponding fragment slides as 
well. 

In explaining the S&S algorithm above, we did not 
prescribe Metropolis sampling for the ends of the total 
path. This does not pose a formal problem for ergodicity, 
but to provide optimal sampling of the path, it is best 
not to unduly neglect any time-slices. In practice, we 
also represent these ends in terms of LC path variables 
which are subjected to Metropolis sampling. 

An important feature of the S&S algorithm is that 
it provides symmetric sampling for all of the time-slices 
in the path, despite the intrinsic asymmetry of the LC 
representation. To illustrate this, consider sampling the 
entire path in the LC representation without fragmenta- 



tion. The time-slice corresponding to the exact middle of 
the path (t/2) is then only moved by changing the a\^\ LC 
path variable. In contrast, the time-slices that directly 
neighbor the middle time-slice will be moved by changing 
LC path variables in any layer. This asymmetry of the 
LC representation causes some path variables to suffer 
from coarse-feature sampling stiffness more drastically 
than others, and it limits the efficiency of sampling the 
entire path to the efficiency of sampling its most coarse 
(and most stiff) layers. The S&S algorithm removes this 
bottleneck by ensuring that a particular time-slice ap- 
pears at all positions within a fragment with equal prob- 
ability. Each time-slice is thus exposed to the full range of 
Monte Carlo moves. Fragmentation, in general, ensures 
multi-resolution sampling of the time-slices, regardless of 
the deficiencies of the underlying random series represen- 
tation. 

Another important feature of the S&S algorithm is that 
it is naturally parallelizable. The Markov property of the 
diffusive dynamics guarantees that the distribution of a 
particular path fragment only depends upon the position 
of its endpoints. By keeping these endpoints fixed during 
the sampling step, we ensure that the fragment distribu- 
tions are statistically independent. It follows that the 
task of sampling the individual fragments can be per- 
formed, in parallel, by up to f+2 (— number of fragments 
+ 2 ends) distributed-memory nodes. We note that since 
a sampling sweep over the path variables in a given frag- 
ment requires nf evaluations of the force field, the cost 
of communication between nodes in the S&S algorithm 
is negligible. 

A consideration in implementing the S&S algorithm 
is the choice of the number of time-slices per fragment. 
In principle, to minimize the total expenditure of com- 
puter hours, choosing nf requires balancing the cost of 
coarse-feature stiffness introduced by the random series 
representation with the cost of (Gibbs-type) correlation 
introduced by excessive fragmentation of the path and 
the corresponding elimination of productive multi-slice 
moves. However, we emphasize that use of the S&S al- 
gorithm with any choice of n,f eliminates the bottleneck 
in sampling efficiency due the asymmetry of the LC rep- 
resentation. It follows that the maximum computational 
penalty to be incurred with excessively long fragments 
is the logarithmic increase in computer time associated 
with attempting inefficient moves in the stiff coarse lay- 
ers. We thus anticipate that the penalty in computa- 
tional cost for overestimating nf in the S&S algorithm is 
a minor concern. This issue is explored in the following 
numerical example. 

In Table [Til we again consider the two-dimensional, 
two-channel numerical example. Using n — 2048 time- 
slices, we sample path distribution using the S&S algo- 
rithm with different numbers of fragments, and we con- 
sider the performance of the algorithm in terms of com- 
puter and wall-clock time. As a measure of the quality 
of sampling in each simulation, we count n^ ops , the num- 
ber of times that the midpoint of the total path crosses 



9 



TABLE II: Sampling efficiency for the two-dimensional exam- 
ple using the parallelized S&S algorithm with various frag- 
ment lengths." 
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a f is the number of fragments, nj + 1 is the number of 
time-slices per fragment, kf is the number of layers in the 
LC series used to describe each fragment. Data in the last 
three columns are normalized. 

barrier in the potential energy surface at x = 0. 

The first three columns of the Table [IT] characterize 
the number and length of the fragments used in the S&S 
algorithm. / is the number of fragments, n/ is the num- 
ber of time-slices per fragment, and kf is the number 
of layers in the LC representation used to describe each 
fragment. In the fourth column, we show the number of 
hopping events observed per number of sweeps through 
the S&S algorithm (each of which includes a Monte Carlo 
sweep through the fragment path variables and a slid- 
ing move). The results are normalized with respect to 
the most coarse fragment. Clearly, better sampling per 
sweep is achieved with larger fragments. This is primar- 
ily due to the fact that increasing the fragment length in 
the LC representation corresponds to attempting addi- 
tional (coarse-layer) Monte Carlo moves per S&S sweep. 
An additional source of improvement is that the correla- 
tion in sampling the fragment endpoints (i.e. Gibbs-type 
correlation) decreases with increasing fragment length. 

Although improved sampling per sweep is obtained 
with longer fragments, coarse-feature stiffness is real- 
ized because these improvements do not fully compen- 
sate for the increasing computational cost of perform- 
ing each sweep. This is illustrated in the fifth column 
of Table HH in which the data from column four are 
scaled by the relative computational costs of perform- 
ing an S&S sweep with different fragment lengths. (To 
compare with the results for the most coarse fragmen- 
tation, each entry in column four is multiplied by the 
factor log 2 (1024)/log 2 (n/) = 10/fc/.) As rt/ increases 
from 16 to 256, we observe the improved sampling ef- 
ficiency that accompanies longer fragments. However, 
at approximately nf = 512, the computational cost of 
performing moves in the coarsest layer exceeds the cor- 
responding improvement in the sampling. Beyond this 
point, we increasingly suffer from the mild coarse- feature 
sampling stiffness of the LC representation. 

The final column in Table |TT] reveals the parallelizabil- 
ity of the S&S algorithm. The number of observed hops 
per unit of wall-clock time is reported by scaling the 



results in column five by the relative number of com- 
putational nodes that can by harnessed via parallcliza- 
tion. (For comparison with the most coarse fragmen- 
tation, each entry in column five is multiplied by the 
factor (/ + 2)/3.) These results again reveal a turnover 
in sampling efficiency - but at much shorter fragmenta- 
tion. Better sampling per unit wall-clock time is gener- 
ally obtained by enlisting the efforts of more computa- 
tional nodes. Eventually, this trend breaks down when 
we lose the benefit of treating the fine-scale features of the 
path with the random series representation. Specifically, 
we see that when we employ fragments with fewer than 
64 time-slices, further parallelization is not beneficial. 

This last observation is consistent with the Ax values 
reported in Table U Both sets of data suggest that, for 
this particular problem and these particular parameters, 
the LC representation is very efficient in sampling path 
variables with fewer than 64 time-slices. Fragmenting the 
path beyond this point diminishes the quality of sam- 
pling. It is convenient that the equilibrated maximum 
displacements indicate the lower bound for the accept- 
able number of time-slices per fragment. 

Before proceeding, we can draw some practical con- 
clusions. The coarse-feature stiffness gives rise to an ex- 
pected turnover in CPU-time efficiency with increasing 
fragment length. However, this penalty is very mild (as 
was also expected), suggesting that little CPU time is 
wasted by utilizing very long fragments. Furthermore, it 
is possible that coarse moves might provide valuable flex- 
ibility in the sampling of complicated problems, regard- 
less of their stiffness. We thus have considerable freedom 
in choosing an acceptable length for the fragments in the 
S&S algorithm. If one only has access to a small num- 
ber of CPU nodes, using the S&S algorithm with long 
fragments promises "optimal" efficiency in terms of com- 
putational effort. On the other hand, if one has access 
to a large number of CPU nodes, then use of shorter 
fragments leads to a lower degree of CPU-time efficiency 
but a vast improvement in wall-clock-time efficiency. To 
summarize: Good sampling efficiency can be expected 
in the S&S algorithm if the fragments are long enough 
to avoid fine-feature sampling stiffness. Beyond this re- 
striction, the final choice of fragment length can be made 
on the basis of the number of available CPU nodes and 
the desire to minimize wall-clock time. Finally, we note 
that the S&S algorithm is sufficiently general to incorpo- 
rate possible improvements in the underlying sampling 
techniques, as they come available. 



V. APPLICATION TO THE 38-ATOM 
LENNARD-JONES CLUSTER 

As a challenging test of our double-ended approach 
to sampling diffusive trajectories, we consider the struc- 
tural rearrangement of the 38-atom Lennard-Jones clus- 
ter (LJ38) at low temperature. Unlike most other small 
Lennard-Jones clusters, LJ38 does not have a global 
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minimum energy structure that exhibits an icosahedral 
corei 22 ' 23 ' 24 ' 25 Instead, its global minimum is a symmet- 
ric truncated octahedron with a total binding energy 
that is only 0.38% lower than the most stable icosahe- 
dral structured The LJ38 potential surface exhibits a 
double-funnel topology, with one basin of stability ex- 
hibiting octahedral-like structures and the other exhibit- 
ing icosahedral-like structures. For temperatures be- 
low T — 0.1 e/ks, the thermal distribution is local- 
ized almost entirely in the octahedral basin. At this 
temperature, a finite-system solid-solid phase transition 
(or "phase change") occurs in which the distribution 
shifts to the entropically favored icosahedral basing At 
T = 0.166 e/kg, a melting transition occurs as the vast 
number of disordered, liquid-like structures become ther- 
mally accessible ! 26 ' 27 

Sampling the LJ38 cluster is known to be extremely 
difficult. 23 Direct Metropolis Monte Carlo sampling of its 
low-temperature thermodynamics is not feasible . 26 ' 28 At 
the low temperatures for which the cluster is solid, ther- 
mally accessible regions of configuration space are sepa- 
rated by very large energetic barriers, and the potential 
surface contains many local minima that pose sampling 
(and kinetic) traps. Furthermore, the octahedral basin, 
and the pathways by which it is accessible, are sufficiently 
narrow to make tortuous any sampling algorithm's task 
of equilibrating the two basins. 

One also expects that the dynamics of interconver- 
sion between the octahedral and icosahedral basins is 
extremely frustrated at low temperature. Our prelimi- 
nary efforts to sample the ensemble of reactive trajec- 
tories using the standard TPS algorithm with shooting 
and shifting confirms this to be the case. Even when long 
trajectories are employed, the system simply can not find 
the basins of stability; it instead becomes trapped in dis- 
ordered, high-energy structures. 

Given the difficulty of ergodically sampling the configu- 
ration space of LJ38 at low temperature, it is unlikely that 
we will be able to fully achieve this feat in path space with 
the techniques introduced here. After all, we have only 
addressed the specific problems associated with sampling 
long dynamical trajectories - not the intrinsic deficiencies 
of the Metropolis Monte Carlo algorithm. However, we 
can expect to harvest reasonable examples of transition 
pathways and to achieve locally-ergodic sampling in the 
space of trajectories. 

We employ the S&S algorithm to sample the distri- 
bution of dynamical trajectories corresponding to inter- 
conversion between the octahedral and lowest icosahe- 
dral structures of the LJ 38 cluster. These trajectories 
connect permutationally distinct isomers of the octahe- 
dral and icosahedral structure. Specifically, we choose 
the two isomers that minimize the Euclidian distance be- 
tween the endpoints of the path at 3.274u, as reported 
by Trygubenko and Wales.— 

The potential energy function, in the Lennard-Jones 



reduced units employed hereafter, is 

38 

V(r) = 4Y / [r^ 12 -rT j 6 }+V c (r), (28) 

i<j 

where is the distance between the i th and j th atoms in 
the cluster, and V c (r) is the constraining potential used 
to prevent evaporation of the cluster. It is defined 

38 

V c (r) = J2v c (5i), (29) 

i=l 

where Si is the distance between the i th atom and the 
center of mass of the cluster, and 

f 0, 5 < r c 

with a = 5000 and r c = 2.65. 

We sample the ensemble of reactive trajectories at the 
phase change temperature, T — 0.1. We employ a fric- 
tion coefficient 7 = 1 and a total simulation time t = 1. 
Two simulations are performed using different fragment 
lengths. In the first, we employ / = 30 fragments of 
length rif = 128 time slices and a number of time-slices 
n = (/+l)*n/ = 3968. In the second, we employ / = 14 
fragments of length rif = 256 and a number of time-slices 
n = 3840. This small difference in n is not sufficient to 
alter the accuracy of sampling. The simulations were 
initialized from a straight line trajectory between the 
minimized icosahedral and octahedral isomers, and these 
endpoints were held fixed during the simulations. The 
simulations were equilibrated for 10 6 S&S sweeps, dur- 
ing which the proposal distribution widths were updated 
to yield approximately 40% acceptance. They were then 
continued for an additional 3 x 10 6 production sweeps. 
Trajectories were recorded every 10 4 sweeps. The sim- 
ulation with n/ = 128 was run in parallel on 32 single- 
processor, distributed-memory nodes, and the simulation 
with rif = 256 employed 16 nodes. The total required 
CPU time, for each simulation, was approximately 10 5 
hours. 

In Table IIIIl we report the equilibrated path variable 
proposal widths for the rif = 256 simulation. The re- 
sults are nearly identical to those obtained for rif = 128, 
and they have the same features as were discussed for 
the two-dimensional example in connection with Table 
UJ In particular, we note a plateau in the Acik.j val- 
ues for k > 6 as the path variables become decoupled 
from the potential energy surface. The plateau value 
seen here is smaller than was found in Table U because 
we are now dealing with three-dimensional, rather than 
two-dimensional, particles. For 3 < k < 6, the de- 
crease with decreasing k while the corresponding Ace val- 
ues remain large, indicating that these layers remain pro- 
ductive in sampling the path. Finally, for k < 2, the Ax 
rapidly decrease as the most coarse layers approach the 
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TABLE III: Proposal distribution widths for the Metropolis 
Monte Carlo moves in the LJ38 application." 
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a k is the layer index for the LC representation of the 
fragment, nt is the number of time-slices moved by changing 
the path variable aj t h, Adfcj is the equilibrated path 
variable distribution width, and Ax is the corresponding 
distance in coordinate space. 



onset of stiffness. This table suggests that a fragmenta- 
tion length of at least nf = 64 is necessary to avoid the 
exclusion of LC path variables that substantially con- 
tribute to the sampling efficiency 

For the purposes of analysis, the sampled trajectories 
were quenched on the potential energy surface using the 
zero temperature string (ZTS) method of E, Ren, and 
Vanden-Eijnden.? Given an initial path, this method 
yields a continuous, minimum-energy path (i.e. one that 
is everywhere tangent to the gradient of the potential 
energy function and is compose of equidistant images of 
the system). We employed the ZTS method using our 
sampled trajectories as input. Specifically, we linearly 
interpolated between every other time-slice in the sam- 
pled diffusive trajectory and performed the ZTS mini- 
mization using a path composed of n/2 images of the 
system. The ZTS calculations were performed in the 
absence of the cluster constraining potential, V c . Con- 
vergence was determined to be reached when the "min- 
max" image, namely the highest energy image along the 
quenched path, exhibited an RMS force of less than 0.001 
per atom. The potential energy profile for a typical tran- 
sition path that has been quenched in this fashion, and 
which illustrates the ruggedness of LJ38 potential energy 
suface, is shown in Figure |4j 

In Fig. [5J we present the time series of minmax val- 
ues obtained from the two S&S simulations at T = 0.1 
reported here. The red and blue lines indicate the re- 
sults for the n/ = 128 and n/ = 256 simulations, re- 
spectively. Trajectories recorded during the equilibration 
stage of the simulation are included. The figure also in- 
cludes two results from a Discrete Path Sampling (DPS) 
calculation^ provided by Wales ~ The solid black line 
indicates the lowest minmax value found to connect the 
permutational isomers, an energy of —168.8346. The 
dashed line indicates the minmax value of —167.7005 
for the most probable path at the lower temperature of 
T = 0.02, according to the harmonic rate theory em- 
ployed in the DPS calculation. 
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FIG. 4: A typical quenched pathway for interconversion be- 
tween the octahedral and icosahedral LJ38 isomers, obtained 
by sampling with the S&S algorithm and quenching using the 
ZTS method. 
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FIG. 5: Minmax values for LJ38 interconversion trajectories 
sampled at T = 0.1 using the S&S algorithm with fragment 
lengths n/ = 128 (red) and rif = 256 (blue). Trajectories are 
output every 10 Monte Carlo sweeps. The solid and dashed 
black lines indicate the minmax values for the lowest energy 
path and most probable path at T=0.02e, respectively, ac- 
cording to a DPS calculation by Wales.— The DPS most 
probable path is determined by a harmonic rate theory.— 



It is clear that the S&S algorithm at least partially 
succeeds in sampling the distribution of paths. It would 
not at all have been surprising if the sampler had sim- 
ply gotten trapped in the region of high-energy path- 
ways around the initial straight-line trajectory. Instead, 
within the first 25 output trajectories, it finds its way 
to lower minmax values and proceeds to explore many 
different transition pathways. It is especially encourag- 
ing that, upon finding a particular low-energy pathway, 
the sampler does not become endlessly trapped. In many 
instances, we see that a low minmax value is found, but 
then abandoned soon after in favor of higher values. This 
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observation also suggests that finite temperature, and 
perhaps the finite trajectory time t, play a role in bi- 
asing the ensemble of trajectories away from the global 
minimum energy path. 

The key role of temperature, and thus entropy, in shap- 
ing the transition path ensemble is corroborated by the 
DPS calculations. Even at the exceedingly low tempera- 
ture of T = 0.02, the lowest energy transition path (solid 
line) is not predicted to be the most favorable (dashed 
line). It is therefore not surprising that our S&S simula- 
tions, which sample trajectory space at the phase change 
temperature of T = 0.1, prefer transition paths that are 
still higher in minmax energy. Both S&S simulations do 
locate the pathway indicated by the dashed line in the 
figure, but they quickly move on in favor of more ener- 
getic paths. These observations highlight the importance 
of methods that actually sample trajectory space. 

Despite these encouraging observations, we can not 
claim that the S&S algorithm achieves fully ergodic sam- 
pling in either simulation. In both S&S simulations re- 
ported in Fig. [5J we see extended periods in which the 
sampler, although not necessarily locally trapped, ap- 
pears to be confined within a "tube" of accessible trajec- 
tories. For example, beyond the first 100 output trajec- 
tories in the n/ = 256 simulation, we see that the system 
consistently returns to a minmax value of —166.675. A 
similar recurrence in the rif — 128 simulation is seen at 
the minmax value of —166.364. These apparent transi- 
tion tubes are analogous to the basins of stability found in 
configuration space, and the sampler's task of equilibrat- 
ing them is difficult for the same reasons. Further char- 
acterization of these transition tubes, perhaps using the 
finite temperature string method in coarse variables^ is 
an interesting possibility for future work. 

Sweep for sweep, the rif = 256 simulation shown in the 
blue should do a better job of sampling than the rif = 128 
simulation shown in the red. The former employs twice 
the fragment length and thus attempts all of the same 
Monte Carlo moves, plus those associated with its most 
coarse LC layer. That the two simulations do not sam- 
ple the same transition tube is an indication that the 
rif = 128, at least, is not ergodic, and most likely neither 
is. However, given the extreme ruggedness of the LJ38 
potential energy surface and the complexity of its low- 
temperature interconversion process, it is very encourag- 
ing to find that reasonable transition paths are obtained 
using the S&S algorithm and that the local sampling of 
path space, at least, is achieved. 

VI. SUMMARY AND CONCLUSIONS 

In this work, we have described several developments 
that aide the efficient transition path sampling of double- 
ended diffusive trajectories. In particular, we have intro- 
duced the symmetric Onsager-Machlup formula, which 



can be utilized as the basis for deriving high-order sym- 
metric integrators that utilize only the forcefield, and we 
have shown that a direct time discretization on a uni- 
form grid generates an approximation that is already ac- 
curate to second order. Furthermore, we have shown that 
the sliding and sampling algorithm can naturally lever- 
age massively parallel computer architectures to simulate 
double-ended trajectories on long timescales. 

Using a two-dimensional example, we demonstrate the 
computational and wall-clock efficiency of the S&S algo- 
rithm. We find that, because of the layered (multi-scale) 
structure of the LC representation and the fact that the 
sliding moves eliminate the intrinsic asymmetries of this 
representation, the computational efficiency of the S&S 
algorithm is relatively insensitive to the length of frag- 
ments employed. This observation greatly simplifies the 
implementation of the algorithm - the fragment lengths 
can be primarily chosen according to the number of avail- 
able nodes and the desire to minimize the wall-clock time 
of simulation. 

We also employ the S&S algorithm to study the low- 
temperature structural interconversion of the LJ38 clus- 
ter. We find that, for the simulation times employed 
here, the S&S algorithm did not ergodically sample the 
full trajectory space - a result that is not at all surprising 
given the difficulty of equilibrating the low-temperature 
thermodynamics of this system with the direct Metropo- 
lis algorithm. However, we are encouraged by finding 
that local sampling of path space is achieved with the 
S&S algorithm and that the sampler was able to find a 
large number of low-energy transition pathways. 

Biological systems at room temperature do not ex- 
hibit the same free energy ruggedness that is found in 
the LJ38 cluster. But they can still have metastable in- 
termediates, entropic bottlenecks, and exceedingly long 
timescales, which frustrate open-ended trajectory sam- 
pling algorithms. The success of the direct Metropo- 
lis Monte Carlo algorithm in sampling the configuration 
space of such systems is a good indication that the S&S 
algorithm will be similarly successful in sampling their 
trajectory space. It is in this regime of biologically rele- 
vant processes that we expect the S&S algorithm to find 
its most immediately applicability. 
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